ending element of the univels is reached where y = y n . Although the integer steps are 
described herein as positive 1, it should be appreciated that the integers can be negative 
and can be in other logical values. It should also be appreciated that the integer values 
that are stored need not correspond to centered values along the primary direction of 
movement. It is just that the centered values provide the most representative sampling 
along the particle track. 

Having stepped to a "next" univel at step 170, it is determined, at step 172, 
whether any of the error terms exceed the threshold value. If the error terms do not 
exceed the threshold values, a determination about the anatomical material of the next 
univel is made at step 174 to see if it is different from the previous or starting univel. 
Again, this is simply done by using the stored integer position values to examine the 
anatomical material mapped in array 152 for that univel against the previous univel. The 
actual points examined are expressed as floats but are only kept track of as integers. 
Thus, as in table 210 (Figure 7), for the next univel having coordinates of (3.83, 2.5, 6.26) 
the anatomical material for that univel is stored in the array at (3,2,6) and a comparison 
between anatomical materials is made against (3,1,5). Similarly, for the univel having 
coordinates (4.13, 3.5, 6.93) the anatomical material for that univel is stored in the array 
at (4,3,6). 

Because of the eventual possibility that stepping in the primary direction of 
movement without stepping along the particle track in the secondary direction of 
movement will cause an error in determining the anatomical material of the univel under 
examination, at step 172, if the error term exceeds the threshold value, an increase in the 
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corresponding coordinate value is performed (step 176) to ensure the proper univel is 
being examined. Thereafter, at step 178, an adjustment of the error terms is performed to 
account for the increase in the corresponding coordinate value. Although not shown, the 
error term could also be adjusted to indicate that stepping only occurred in the primary 
direction of movement. Thence, once adjusted, the determination of the anatomical 
material of the next univel is made at step 174. 

It should be appreciated that the anatomical material of the "next" univel is made 
in comparison to the starting univel, or, as the movement of the particle is tracked along 
the particle track, is made in comparison to the previous univel. If, at step 174, the 
anatomical material is not at least substantially different, the movement of the particle 
along the particle track is reiteratively traversed to the next univel (step 170) until 
eventually the particle exits the geometry or intersects with a new material. 

Thus, at step 174, if the anatomical material of the next univel is different from 
the previous or starting univel, a determination is made, at step 180, to see if the particle 
has exited the geometric model. As in the prior art, if the particle has exited the 
geometric model, the particle transport simulation is terminated at step 182. 

If the particle has not exited the geometric model at step 180 and the anatomical 
material of the next univel is different from the previous univel, the position of 
intersection with the new material is determined at step 184. When determined, this 
position of intersection is reported at step 186 for use in another part of the computer 
executable instructions. 
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It should be appreciated from the foregoing that computational time is greatly 
preserved by stepping through the geometric model in integer based increments because 
each of the stepping computations and each determination about the anatomical material 
of each univel is performed by the computer without requiring the use of floating point 
mathematics. Thus, a medical image having pixels of information in 512 x 512 
resolution x 5 12 axial slices, millions of computations are performed over the course of 
numerous particles emanated from a radiation source. As described subsequently, this 
reduction in tracking time has been shown to be at least one order of magnitude faster 
than any computations heretofore known in the field. 

The step 184 for determining where the position of intersection with the new 
material happens, is further described with reference to Figures 8 and 9. In Figure 8, it is 
known that in some univel 148, the particle traveling along the particle track entered an 
anatomical material different from the previous univel. To determine the precise 
intersection, it is first known that the particle entered the univel 148 along the particle 
track through one of three planes. The particle may have entered the univel: through the 
X-Y plane as along particle track 220; through the X-Z plane along particle track 222; or 
through the Y-Z plane along particle track 224. The X, Y and Z planes being taken in 
reference to the Cartesian coordinate system depicted. Again, other coordinate systems 
can be used. 

In a preferred embodiment, with reference to Figure 9, three possible intersection 
points are established along the three primary planes described above, step 190. For 
example, the first position is 221 along particle track 220. The second position is 223 

DocketNo.LIT-PI-316 

29 



along particle track 222. The third position is 225 along particle track 224. Since each of 
these three positions are along a planar surface of a univel, a small epsilon may be added 
to move each of the three positions inside the univel by a small amount so that ambiguity 
of being on the planar surface can be avoided for computational purposes. 

In a preferred embodiment, the three positions correspond to a floor or ceiling 
operator. The floor or ceiling is in reference to whether the particle track is moving in 
positive or negative increments. If positive, a floor is set. If negative, a ceiling is set. An 
example of this is depicted by particle track 222 in which positive increment 
advancement occurs in the Y and Z directions and negative increment advancement 
occurs in the X direction. Then, at step 192, the position where the particle track first 
enters the new material is determined. This is done by examining whichever particle 
through one of the three positions first hits or intersects the anatomical material that is 
different from the previous univel, this is the position of intersection. Again, this 
intersection position is reported at step 186. Floor and ceiling operators are well known 
in the art and are not described herein in detail. 

With this method of integer based tracking of a movement of a particle along a 
particle track through a geometric model, it should be appreciated that some univels may 
be skipped over when tracking the particle. An example of this is shown with reference 
to Figure 10, wherein a particle track 230, shown only in the X-Y plane, traverses through 
a small corner of univel 232. As such, if univel 232 is of the same anatomical material as 
univel 234, there is no need to perform a detailed examination regarding this univel and 
progression of the particle track can continue to univel 236. Thus, it is only when a 
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univel has a different anatomical material from the previous univel that any further 
detailed calculations are required to be made. If the anatomical material of univel 232 is 
different, but the particle track reenters the original material in univel 236, then an 
insufficient volume in 234 was intersected to count as a boundary crossing. In the case 
where the anatomical material of 236 is different than 234, univel 232 will also be 
examined when determining the precise crossing into the new material since three planes 
of entry into the new material are considered. Again, when calculations for particle 
tracks are performed through millions of univels, tremendous computational time is 
saved. 

In contrast to the prior art, it should be appreciated that computational accuracy is 
improved with more representations of the treatment volume than with fewer. For 
example, some methods in the prior art used 500 pixels as a single element representation 
for tracking a particle movement. Yet, in a 256 x 256 resolution, this only equates to 
about 130 elements in the model. If a particular particle track only passed through a 
small portion of these 130 elements, an accurate understanding required for 
computational dosimetry would be severely lacking. Yet here, a 256 x 256 resolution 
equates to 65,536 univels per axial slice. But because the tracking is performed in 
integer based increments, the tracking is not only faster but yields more accurate data in 
the dosimetry planning. 

In an alternate embodiment, after step 166 (Figure 6 A), a decision 167 is made 
whether to follow along the particle track or not. When performing Monte Carlo 
simulation using an alternative scheme known as "boundary elimination," it is only 
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necessary to know the material of the starting univel and not required to follow along the 
particle track to determine the next material intersection. Thus, for this alternate method, 
and for some editing purposes, return is made to the calling program immediately after 
determining the material of the starting univel. As such, this alternative step is indicated 
by dashed lines. 

With reference to Figure 1 1, it should be appreciated that as medical imagery 
becomes even more sophisticated, it is expected that even greater resolutions will be 
provided, such as in a 1000 pixel x 1000 pixel resolution with 1000 axial slices. Thus, to 
improve computational times for tracking a particle movement through the geometric 
model, groupings of elements may be advantageously arranged. One such grouping uses 
a super univel 250 comprised of an arrangement of smaller univels in a 2 x 2 x 2 
configuration. Still other combinations of univels can be effectuated. 

Example 1 

The following represents data obtained from tracking a movement of 100,000 
particles along random particle tracks (Monte Carlo simulated) through a geometric 
model constructed from a 256 x 256 x 33 medical image consisting of a buffer material, 
scalp, skull, brain and tumor anatomical materials. 

The particle tracks began at a random initial position in the geometric model and 
traversed in a random direction. Each movement was tracked along the particle track 
until either the particle intersected an anatomical material different from the anatomical 
material of the previous univel or was exited from the geometric model. Of the 100,000 
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particle tracks, 55,137 positions of intersections and 44,863 exits from the geometric 
model were reported. 



Table 1 



3 600 422 

, \J\J V)"^** 


univels having particle tracks 


33,670.034 


positions of intersection/ sec 


1,212,263.3 


univels/ sec 


36.004 


univels tracked through/ position of 
intersection reported 


2.970 


elapsed time (sec) 


196,270.562 


distance traveled all particle tracks (cm) 


66,084.364 


distance traveled/ sec 



It should be appreciated that since the simulated particle transport was performed 
in less than about 0.2 hours, the foregoing represents an advance over the present state of 
the art by as much as 51 times. Heretofore, such simulated particle transport would 
routinely require as much as 10 hours of computational time or more. 

Example 2 

The following represents the actual algorithm information used to simulate such 
advanced particle transport along a particle track as presented at the 1998 Radiation 
Protection and Shielding Division Topical Conference in Nashville, TN in April. Note 
that the subject matter of this presentation was directed only to external neutron sources. 

Data Initialization: The uniformly spaced medical image data is read into an 
array. The x-pixel-size, y-pixel-size, and z-pixel size along with the minimum value of 
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each coordinate is stored so conversions between world coordinates (WC) and normalized 
array coordinates (NAC) can be easily made. Here, the NAC simply corresponds to a 
location in the array of univels. For example, any location in the array can be found by 
an ordered triple of nonnegative integers, i.e., lookup (x,y,z) = array (z (width x length) + 
y (width) + x). A univel in WC is A mm x B mm x C mm. Whereas the univel in NAC 
is 1 x 1 x 1, for example. 

Parameters: A call to the movement of the particle along a particle track is of 

the form: 

Track_Ray (position_vector, direction_unit_vector, ptrJo_miss_flag, 



ptr_to_current_region, ptr_to_next_region, 



ptr_to_distance_to_next_region); 



Input to algorithm: 



position_vector: 



Initial position of particle track in WC 



direction unit vector: Normalized direction of particle track in WC 



Output of algorithm: 



miss_flag: 



Either hit a new region or exit the geometric 



model 



current_region: 



The region (univel) the particle track starts 



m 



next_region: 



The first region intersected 



distance_to_next_region: 



The distance to the next_region (univel) in 



WC 
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Algorithm Initialization Calculations: The initial position and direction must 
be converted from WC to NAC. The initial anatomical material is stored in 
current_region. If the particle track does not start inside the univel geometric model, an 
intersection point with the univel geometric model must be calculated and an artificial 
starting point is set at this boundary intersection with the outer univel. 

Stepping Algorithm: Though the internal routines of the algorithms vary, each is 
based on using integer arithmetic to find univels that the ideal particle track passes 
through. Each investigated univel has a corresponding call to a function that looks up the 
anatomical material type of the univel at the given position. The stepping algorithm 
terminates when a univel of a new anatomical material type is found or the particle along 
the particle track exits the geometric model. 

Algorithm Completion Calculations: The position of intersection is computed 
accurately or miss_flag is set to indicate the particle exited the geometric model without 
an intersection. The distance to this point is calculated in WC and returned in 
distance_to_next_region. The new material encountered is stored in next_region. 

Example 3 

The following data was presented at the 1998 conference in Nashville, TN and is 
exemplary of a particle track having Y as a primary direction of movement, X is the 
secondary direction of movement increasing in 0.125 units of a Cartesian coordinate 
system and Z is the tertiary direction of movement increasing in 0.75 units. The initial 
starting position of the particle track after centering is x 0 = 5.00, y 0 = 1.5 and z 0 = 10.125. 
Truncating (trunc) is the rounding down function. Stepping along the primary direction 
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of movement yields the following data with an error term being an integer in the range of 
-32,768 to 32,767: 

Table 2 



A 


v 

y 


z 


trunc(x) 


trunc(y) 


trunc(z) 


J.UUU 




10.125 


5 


1 


10 


J. 1ZJ 




10.875 


5 


2 


10 


5.250 


3.5 


11.625 


5 


3 


11 


5.375 


4.5 


12.375 


5 


4 


12 


5.500 


5.5 


13.125 


5 


5 


13 


5.625 


6.5 


13.875 


5 


6 


13 


5.750 


7.5 


14.625 


5 


7 


14 


5.875 


8.5 


15.375 . 


5 


8 


15 


6.000 


9.5 


16.125 


6 


9 


16 



The bulk of the corresponding stepping algorithm for this example is as follows: 
ADDX = 0.125*32768 = 4096 
ADDZ = 0.750*32768 = 24576 

ADDXJDECERR = ADDX - 32768 
ADDZDECERR = ADDZ - 32768 

ERRX = (x 0 - trunc(x 0 ))*32768 + ADDXJDECERR = -28672 
ERRZ = (z 0 - trunc(z 0 ))*32768 + ADDZ DECERR = -4096 

XI = trunc(x) 
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YI = trunc(y) 
ZI = trunc(z) 

BEGIN_LOOP 
LOOKUP(XI,YI,ZI) 
YI = YI+1 
If(ERRX>=0) 
XI = XI+1 

ERRX = ERRX + ADDXDECERR 
Else 

ERRX = ERRX + ADDX 
If(ERRZ>=0) 
ZI = ZI + 1 

ERRZ = ERRZ + ADDZDECERR 
Else 

ERRZ = ERRZ + ADDZ 
END_LOOP 

The next table shows the values computed by the algorithm. Notice that the error 
term, i.e., ERRX or ERRZ, is a pre-computation used to determine how XI and ZI will 
change in the next iteration, increasing by 1 if the error is greater than or equal to 0 and 
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remaining the same otherwise. The steps are similar when the directions are allowed to 
be decreasing. 

Table 3 



VT 
Al 


YT 

X 1 


ZI 


ERRX 


ERRZ 


J 


1 

1 


10 


-28672 


-4096 


c 


o 


10 


-24576 


20480 


5 


3 


11 


-20480 


12288 


5 


4 


12 


-16384 


4096 


5 


5 


13 


-12288 


-4096 


5 


6 


13 


-8192 


20480 


5 


7 


14 


-4096 


12288 


5 


8 


15 


0 


4096 


6 


9 


16 


-28672 


-4096 



This example was for each x, y and z increasing and y varying the most. The 
method would be similar for either x or z varying the most. Negative directions only 
cause minor complications wherein absolute values of the directions are used to compute 
the ADDx's. ERRx's are computed based on the distance to the ceiling side rather than 
the floor side of the initial univel from the initial point. 

In the general case, there is some roundoff error since the initial positions and 
increments may not be expressed exactly as fractions of 32768. Using 32 bit arithmetic 
instead, an individual error is less than 2" 31 (=4.66* 10" 10 ) and the error accumulates by at 
most that much on each iteration. In a 256x256x40 array of univels, the cumulative error 
could at worst be 256*2" 31 = 2" 23 (* 1 .19*10" 7 ) which is insignificant in most cases for two 
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reasons. The algorithm yields only an approximate x, y, and z as an array position which 
is then refined to give a precise intersection that is not subject to this cumulative error. 
Also, the approximated particle movement follows very closely the same univels as the 
ideal particle track. Being off by at most 2" 23 of a side length suggests that the 
approximate particle track reports a position of intersection not intersected by the ideal 
particle track less than 1 in 1,000,000 times. If the algorithm reports a position of 
intersection that is not verifiable, the particle movement along the particle track is simply 
allowed to continue. Any position of intersection distance returned is precise and 
verified. 

The present invention may be embodied in other specific forms without departing 
from its spirit or essential characteristics. The described embodiments are to be 
considered in all respects only as illustrative and not restrictive. The scope of the 
invention is, therefore, indicated by the appended claims rather than by the foregoing 
description. All changes which come within the meaning and range of equivalency of the 
claims are to be embraced within their scope. 
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